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Abstract. This work deals with the position control of selected patterns in reaction- 
diffusion systems. Exemplarily, the Schlogl and FitzHugh-Nagumo model are dis¬ 
cussed using three different approaches. First, an analytical solution is proposed. 
Second, the standard optimal control procedure is applied. The third approach ex¬ 
tends standard optimal control to so-called sparse optimal control that results in 
very localized control signals and allows the analysis of second order optimality 
conditions. 


Introduction 

Beside the well-known Turing patterns, reaction-diffusion (RD) systems pos¬ 
sess a rich variety of self-organized spatio-temporal wave patterns including 
propagating fronts, solitary excitation pulses, and periodic pulse trains in one¬ 
dimensional media. These patterns are “building blocks” of wave patterns like 
target patterns, wave segments, and spiral waves in two as well as scroll waves 
in three spatial dimensions, respectively. Another important class of RD pat¬ 
terns are stationary, breathing, and moving localized spots [1-7]. 

Several control strategies have been developed for purposeful manipula¬ 
tion of wave dynamics as the application of closed-loop or feedback-mediated 
control loops with and without delays [8-11] and open-loop control that in¬ 
cludes external spatio-temporal forcing [10,12-14], optimal control [15-17], 
and control by imposed geometric constraints and heterogeneities on the 
medium [18,19]. While feedback-mediated control relies on continuously moni¬ 
toring of the system’s state, open-loop control is based on a detailed knowledge 
of the system’s dynamics and its parameters. 

Experimentally, feedback control loop have been developed for the photo¬ 
sensitive Belousov-Zhabotinsky (BZ) reaction. The feedback signals are ob¬ 
tained from wave activity measured at one or several detector points, along 
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detector lines, or in a spatially extended control domain including global feed¬ 
back control [8,9,20]. Varying the excitability of the light-sensitive BZ medium 
by changing the globally applied light intensity forces a spiral wave tip to de¬ 
scribe a wide range of hypocycloidal and epicycloidal trajectories [21,22]. 
Moreover, feedback-mediated control loops have been applied successfully in 
order to stabilize unstable patterns in experiments, such as unstable travel¬ 
ing wave segments and spots [11]. Two feedback loops were used to guide 
unstable wave segments in the BZ reaction along pre-given trajectories [23]. 
An open loop control was successfully deployed in dragging traveling chemi¬ 
cal pulses of adsorbed CO during heterogeneous catalysis on platinum single 
crystal surfaces [24]. In these experiments, the pulse velocity was controlled 
by a laser beam creating a movable localized temperature heterogeneity on 
an addressable catalyst surface, resulting in a V-shaped pattern [25]. Drag¬ 
ging a one-dimensional chemical front or phase interface to a new position by 
anchoring it to a movable parameter heterogeneity, was studied theoretically 
in [26,27]. 

Recently, an open-loop control for controlling the position of traveling 
waves over time according to a prescribed protocol of motion cj){t) was pro¬ 
posed that preserves simultaneously the wave profile [28]. Although position 
control is realized by external spatio-temporal forcing, i.e., it is an open- 
loop control, no detailed knowledge about the reaction dynamics as well as 
the system parameters is needed. We have already demonstrated the ability 
of position control to accelerate or decelerate traveling fronts and pulses in 
one spatial dimension for a variety of RD models [29,30]. In particular, we 
found that the analytically derived control function is close to a numerically 
obtained optimal control solution. A similar approach allows to control the 
two-dimensional shape of traveling wave solutions. Control signals that real¬ 
ize a desired wave shape are determined analytically from nonlinear evolution 
equations for isoconcentration lines as the perturbed nonlinear phase diffu¬ 
sion equation or the perturbed linear eikonal equation [31]. In the work at 
hand, we compare our analytic approach for position control with optimal 
trajectory tracking of RD patterns in more detail. In particular, we quantify 
the difference between an analytical solution and a numerically obtained re¬ 
sult to optimal control. Thereby, we determine the conditions under which 
the numerical result approaches the analytical result. This establishes a basis 
for using analytical solutions to speed up numerical computations of optimal 
control and serves as a consistency check for numerical algorithms. 

We consider the following controlled RD system 

dtu{x, t) — 'DAu{x, t) + R{u{x, t)) = Bf{x, t). (1) 

Here, u{x^t) = (i4i(x,t),... ,Un{x,t))'^ is a vector of n state components in 
a bounded or unbounded spatial domain i? C of dimension N G {1, 2, 3}. 
V is an n X n matrix of diffusion coefficients which is assumed to be diagonal, 
X> = diag(Di,... ,Dn), because the medium is presumed to be isotropic. A 
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represents the A^-dimensional Laplacian operator, and R denotes the vector 
of n reaction kinetics which, in general, are nonlinear functions of the state. 
The vector of control signals /(x,t) = (/i(x, t),...,/^(x, t))^ acts at all 
times and everywhere within the spatial domain Q. The latter assumption 
is rarely justified in experiments, where the application of control signals is 
often restricted to subsets of Q. However, notable exceptions, as e.g. the al¬ 
ready mentioned photosensitive BZ reaction, exist. Here, the light intensity is 
deployed as the control signal such that the control acts everywhere within a 
two-dimensional domain. 

Equation (1) must be supplemented with an initial condition ii(x,to) = 
Uq{x) and appropriate boundary conditions. A common choice are no-flux 
boundary conditions at the boundary U = dQ x (0,T), dnu{x^t) = 0 , where 
dnU denotes the component-wise spatial derivative in the direction normal to 
the boundary F = dQ of the spatial domain. 

Typically, the number m of independent control signals in Eq. (1) is smaller 
than the number n of state components. We call such a system an underac¬ 
tuated system. The n x m matrix B determines which state components are 
directly affected by the control signals. If m = n and the matrix B is regular, 
it is called a fully actuated system. 

Our main goal is to identify a control / such that the state u follows a 
desired spatio-temporal trajectory Ud^ also called a desired distribution, as 
closely as possible everywhere in space Q and for all times 0 < t < T. We 
can measure the distance between the actual solution u of the controlled RD 
system Eq. (1) and the desired trajectory Ud up to the terminal time T with 
the non-negative functional 

J{u) =\\u-Ud\\l2^Q), ( 2 ) 

where || • 11^2 is the L^-norm defined by 

T 

whUq) = JJ dxdt (3) 

in the space-time-cylinder Q := Q x (0,T). The functional Eq. (2) reaches its 
smallest possible value, J = 0, if and only if the controlled state u equals the 
desired trajectory almost everywhere in time and space. 

In many cases, the desired trajectory Ud cannot be realized exactly by the 
control, cf. Ref. [32] for examples. However, one might be able to find a control 
which enforces the state u to follow Ud as closely as possible as measured by 
J. A control / = / is optimal if it realizes a state u which minimizes J. The 
method of optimal control views J as a constrained functional subject to u 
satisfying the controlled RD system Eq. (1). 

Often, the minimum of the objective functional J, Eq. (2), does not exist 
within appropriate function spaces. Consider, for example, the assumption 
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that the controlled state, obtained as a solution to the optimization problem, 
is continuous in time and space. Despite that a discontinuous state u leading 
to a smaller value for J{u) than any continuous function might exist, this 
state is not regarded as a solution to the optimization problem. Furthermore, 
a control enforcing a discontinuous state may diverge at exactly the points of 
discontinuity; examples in the context of dynamical systems are discussed in 
Ref. [32]. For that reason, the unregularized optimization problem, Eq. (2), 
is also called a singular optimal control problem.To ensure the existence of a 
minimum of J and bounded control signals, additional (inequality) constraints 
such as bounds for the control signal can be introduced, cf. Ref. [33]. Alter¬ 
natively, it is possible to add a so-called Tikhonov-regularization term to the 
functional Eq. (2) which is quadratic in the control, 

J(U, /) = \\U - 'Ud|||2(Q) + T^\\f\\l2(^Qy (4) 

The I/^-norm of the control / is weighted by a small coefficient z/ > 0. This 
term might be interpreted as a control cost to achieve a certain state u. 
Since the control / does not come for free, there is a “price” to pay. In 
numerical computations, z/ > 0 serves as a regularization parameter that 
stabilizes the algorithm. Eor the numerical results shown in later sections, we 
typically choose z/ in the range 10“^ < < 10“^. While z/ > 0 guarantees 

the existence of an optimal control / in one and two spatial dimensions even 
in the absence of bounds on the control signal [33], it is not known whether 
Tikhonov-regularization alone also works in spatial dimensions N larger than 
two. Here, we restrict our investigations to one and two spatial dimensions. 
The presence of the regularization term causes the states to be further away 
from the desired trajectories than in the case of z/ = 0. Thus, the case z/ = 0 
is of special interest. Naturally, the solution n for z/ = 0 is the closest (in 
the L‘^ {Q)-sense) to the desired trajectory Ud among all optimal solutions 
associated with any z/ > 0. Therefore, it can be seen as the limit of realizability 
of a certain desired trajectory Ud. 

In addition to the weighted L^-norm of the control, other terms can be 
added to the functional Eq. (4). An interesting choice is the weighted L^-norm 
such that the functional reads 

j(m, /) = i|m - Mdiii2(Q) + 2^ii/iii2(Q) + k||/i|li(q)- (5) 

Eor appropriate values of />: > 0, the corresponding optimal control becomes 
sparse, i.e., it only acts in some localized regions of the space-time-cylinder 
while it vanishes identically everywhere else. Therefore, it is also called sparse 
eontrol or sparse optimal eontrol in the literature, see Refs. [34-37]. In some 
sense, we can interpret the areas with non-vanishing sparse optimal control 
signals as the most sensitive areas of the RD patterns with respect to the 
desired control goal. A manipulation of the RD pattern in these areas is most 
efficient while control signals applied in other regions have only weak impact. 
Eurthermore, the weighted L^-norm enables the analysis of solutions with a 



Analytical and Optimal Position Control of Reaction-Diffusion patterns 


5 


Tikhonov-regularization parameter v tending to zero. This allows to draw 
conclusions about the approximation of solutions to unregularized problems 
by regularized ones. 

In Sect. 1, we present an analytical approach for the control of the posi¬ 
tion of RD patterns in fully actuated systems. These analytical expressions 
are solutions to the unregularized {y = 0) optimization problem, Eq. (2), 
and might provide an appropriate initial guess for numerical optimal con¬ 
trol algorithms. Notably, neither the controlled state nor the control signal 
suffering from the problems are usually associated with unregularized opti¬ 
mal control; both expressions yield continuous and bounded solutions under 
certain assumptions postulated in Sect. 1. In Sect. 2, we state explicitly the 
optimal control problem for traveling wave solutions to the Schlogl [1,38] 
and the FitzHugh-Nagumo model [39,40]. Both are well-known models to de¬ 
scribe traveling fronts and pulses in one spatial dimension, solitary spots and 
spiral waves in two spatial dimensions, and scroll waves in three spatial di¬ 
mensions [4,10,41,42]. We compare the analytical solutions from Sect. 1 with 
a numerically obtained regularized optimal control solution for the position 
control of a traveling front solution in the one-dimensional Schlogl model in 
Sect. 2.3. In particular, we demonstrate the convergence of the numerical re¬ 
sult to the analytical solution for decreasing values v. The agreement becomes 
perfect within numerical accuracy if v is chosen sufficiently small. Section 3 
discusses sparse optimal control in detail and presents numerical examples 
obtained for the FitzHugh-Nagumo system. Finally, we conclude our findings 
in Sect. 4. 


1 Analytical approach 

Below, we sketch the idea of analytical position control of RD patterns 
proposed previously in Refs. [28,29]. For simplicity, we consider a single¬ 
component RD system of the form 

dtu{x, t) — Dd^u{x^ t) + R{u{x, t)) = /(x, t), (6) 

in a one-dimensional infinitely extended spatial domain x G M. The state u 
as well as the control signal / are scalar functions and the system Eq. (6) 
is fully actuated. Usually, Eq. (6) is viewed as a differential equation for the 
state u with the control signal / acting as an inhomogeneity. Alternatively, 
Eq. (6) can also be seen as an expression for the control signal. Exploiting 
this relation, one simply inserts the desired trajectory for u in Eq. (6) and 
obtains for the control 

f{x, t) = dtUd{x, t) - Dd^Udix, t) + R{ud{x, t)). (7) 

In the following, we assume that the desired trajectory Ud is sufficiently 
smooth everywhere in the space-time-cylinder Q such that the evaluation of 
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the derivatives dtUd and d‘^Ud yields continuous expressions. We call a desired 
trajectory Ud exactly realizable if the controlled state u equals Ud everywhere in 
Q, i.e., u{x^ t) = Ud{x^ t). For the control signal given by Eq. (7), this can only 
be true if two more conditions are satisfied. First, the initial condition for the 
controlled state must coincide with the initial state of the desired trajectory, 
i.e., u{x,to) = Ud{x^to). Second, all boundary conditions obeyed by u have to 
be obeyed by the desired trajectory Ud as well. Because of u{x,t) = Ud{x,t), 
the corresponding unregularized functional J, Eq. (2), vanishes identically. 
Thus, the control / is certainly a control which minimizes the unregularized 
functional J and, in particular, it is optimal. 

In conclusion, we found a solution to the unregularized optimization prob¬ 
lem Eq. (2). The solution for the controlled state is simply u{x,t) = Ud{x,t), 
while the solution for the control signal is given by Eq. (7). Even though we 
are dealing with an unregularized optimization problem, the control signal as 
well as the controlled state are continuous and bounded functions provided 
the desired trajectory Ud is sufficiently smooth in space and time. 

Generalizing the procedure to multi-component RD systems in multiple 
spatial dimensions, the expression for the control reads 

-'DAud{x,t) + R{ud{x,t))). (8) 

Once more, the initial and boundary conditions for the desired trajectory 
Ud have to comply with the initial and boundary conditions of the state u. 
Clearly, the inverse of B exists if and only if B is a regular square matrix, i.e., 
the system must be fully actuated. We emphasize the generality of the result. 
Apart from mild conditions on the smoothness of the desired distributions Ud, 
Eq. (8) yields a simple expression for the control signal for arbitrary Ud- 

Next, we consider exemplarily the position control of traveling waves (TW) 
in one spatial dimension. Traveling waves are solutions to the uncontrolled 
RD system, i.e., Eq. (1) with / = 0. They are characterized by a wave profile 
u{x,t) = Uc{x — ct) which is stationary in a frame of reference ^ = x — ct 
co-moving with velocity c. The wave profile Uc satisfies the following ordinary 
differential equation (ODE), 

D^7"(C) + c^7'(C)-i^(^7e(e)=0, (9) 

The prime denotes differentiation with respect to Note that stationary 
solutions with a vanishing propagation velocity c = 0 are also considered as 
traveling waves. The ODE for the wave profile, Eq. (9), can exhibit one or more 
homogeneous steady states. Typically, the wave profile Uc approaches either 
two different steady states or the same steady state for ^ ^ Too. This fact 
can be used to classify traveling wave profiles. Front profiles connect different 
steady states for ^ ^ Too and are found to be heteroclinic orbits of Eq. (9). 
Pulse profiles join the same steady state and are found to be homoclinic 
orbits [43]. Furthermore, all TW solutions are localized in the sense that their 
spatial derivatives of any order m > 1 decay to zero, lim^^ioo d^Uc{^) = 0. 
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We propose a spatio-temporal control signal which shifts the trav¬ 

eling wave according to a prescribed protocol of motion 0(t) while simul¬ 
taneously preserving the uncontrolled wave profile t/c- Correspondingly, the 
desired trajectory reads 

Ud{x, t) = Uc{x - 0(t)). (10) 

Note that the desired trajectory is localized for all values of 0(t) because the 
TW profile Uc is localized. The initial condition for the state is u{x,to) = 
Uc{x — xo) which fixes the initial value of the protocol of motion as 4>{to) = xq. 
Then, the solution Eq. (8) for the control signal becomes 

f{x,t) = B-\-4>{t)U'dx - - VU'^ix - cPit)) + R{Udx - 4>m), (11) 

with denoting the derivative of 0(t) with respect to time t. Using Eq. (9) 
to eliminate the non-linear reaction kinetics R, we finally obtain the following 
analytical expression for the control signal 

fix, t) = ic- m)B-^U',ix - m) =■ /an. (12) 

Remarkably, any reference to the reaction function R drops out from the ex¬ 
pression for the control. This is of great advantage for applications without or 
only incomplete knowledge of the underlying reaction kinetics R. The method 
is applicable as long as the propagation velocity c is known and the uncon¬ 
trolled wave profile Uc can be measured with sufficient accuracy to calculate 
the derivative U'^. 

Being an open loop control, a general problem of the proposed position 
control is its possible inherent instability against perturbations of the initial 
conditions as well as other data uncertainty. However, assuming protocol ve¬ 
locities (j){t) close to the uncontrolled velocity (j) ^ the control signal 
Eq. (12) is small in amplitude and enforces a wave profile which is relatively 
close to the uncontrolled TW profile Uc- Since the uncontrolled TW is pre¬ 
sumed to be stable, the controlled TW might benefit from that and a stable 
open loop control is expected. This expectation is confirmed numerically for 
a variety of controlled RD systems [28], and also analytically in Ref. [29]. 

Despite the advantages of our analytical solution stated above, there are limits 
for it as well. The restriction to fully actuated systems, i.e., systems for which 
exists, is not always practical. In experiments with RD systems, the 
number of state components is usually much larger than one while the number 
of control signals is often restricted to one or two. Thus, the question arises 
if the approach can be extended to underactuated systems with a number 
of independent control signals smaller than the number of state components. 
This is indeed the case but entails additional assumptions about the desired 
trajectory. In the context of position control of TWs, it leads to a control 
which is not able to preserve the TW profile for all state components, see 
Ref. [28]. The general case is discussed in the thesis [32] and is not part of 
this paper. 
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Moreover, in applications it is often necessary to impose inequality con¬ 
straints in form of upper and lower bounds on the control. For example, the 
intensity of a heat source deployed as control is bounded by technical reasons. 
Even worse, if the control is the temperature itself it is impossible to attain 
negative values. Since the control signal /an foi* position control is propor¬ 
tional to the slope of the controlled wave profile the magnitude of the 
applied control may locally attain non-realizable values. In our analytic ap¬ 
proach no bounds for the control signal are imposed. The control signal / as 
given by Eq. (8) is optimal only in case of a vanishing Tikhonov-regularization 
parameter z/ = 0, cf. Eq. (4). Moreover, desired trajectories Ud which do not 
comply with initial as well as boundary conditions or are non-smooth might be 
requested. Lastly, the control signal / cannot be used in systems where only 
a restricted region of the spatial domain Q is accessible by control. While all 
these cases cannot be treated within the analytical approach proposed here, 
optimal control can deal with many of these complications. 


2 Optimal Control 

In the following, we recall the optimal control problem and sketch the most 
important analytical results to provide the optimality system. 

2.1 The Control Problem 

Eor simplicity, we state the optimal control problem explicitly for the EitzHugh- 
Nagumo system [39,40]. The EitzHugh-Nagumo system is a two-component 
model u = for an activator u and an inhibitor 

dtu{x, t) — Au{x, t) -h R{u{x, t)) -h a v{x, t) = /(x, t), 
dtv{x, t) (3 v{x, t) — ^u{x, t) S = 0, 

in a bounded Lipschitz-domain i? C of dimension I < N < 
the single-component control / appears solely on the right-hand side of the 
first equation, this system is underactuated. Allowing a control in the second 
equation is fairly analogous. The kinetic parameters a, 7 , and 6 are real 
numbers with P > 0. Moreover, the reaction kinetics is given by the nonlinear 
function R{u) = u(u — a)(u — 1) for 0 < a < 1. Note that the equation for 
the activator u decouples from the equation for the inhibitor -z; for a = 0 , 
cf. Eqs. (13), resulting in the Schlogl model [1,38], sometimes also called the 
Nagumo model. We assume homogeneous Neumann-boundary conditions for 
the activator u and iz(x, 0 ) = uq{x)^ v{x^Q) = V{){x) are given initial states 
belonging to i.e., they are bounded. 

The aim of our control problem is the tracking of desired trajectories 
Ud = {ud^VdY" in the space-time cylinder Q and to reach desired terminal 
states Ut = at the final time T. In contrast to the analytic ap¬ 

proach from Sect. 1 , these desired trajectories are neither assumed to be 


(13) 
3. Since 
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smooth nor compatible with the given initial data or boundary conditions. 
For simplicity, we assume their boundedness, i.e., G and 

{ut^vt)'^ G The goal of reaching the desired states is expressed as 

the minimization of the objective functional 

J(u,v,f) = 1 [ct\\u{ - ,T) - Ut\\12(q) + cX.\\v{ - ,T) - 

f (14) 

+ 2 (^d l|w - «d|lL2(Q) + Cd Ik - ^'d|lL2(Q)) + ^II/IIl2(Q)- 

This functional is slightly more general than the one given by Eq. (2) because 
it also takes into account the terminal states. We emphasize that the given 
non-negative coefficients c^, c^, c^, and can also be chosen as functions 
depending on space and time. In some applications, this turns out to be very 
useful [44]. The control signals can be taken out of the set of admissible 
controls 


J^a.d = {f & L°°{Q) ■ fa < f{x,t) < h, ioi{x,t)eQ}. ( 15 ) 

The bounds —oo<fa<fb<oo model the technical capacities for generating 
controls. 

Under the previous assumptions, the controlled RD equations (13) have a 
unique weak solution denoted by (uf^Vf)'^ for a given control / G J^ad- This 
solution is bounded, i.e., Uf, Vf G cf. [44]. If the initial data 

are continuous then Uf and Vf are continuous on i? x [0, T] with i? = i? U df2 
as well. Moreover, the control-to-state mapping G := / i-G {uf^Vf)'^ is twice 
continuously (Frechet-) differentiable. A proof can be found in Ref. [44, The¬ 
orem 2.1,Corollary 2.1, and Theorem 2.2]. Expressed in terms of the so¬ 
lution {uf^Vf)'^, the value of the objective functional depends only on /, 
J(u^v^f) = J{uf,Vf,f) =: F{f)^ and the optimal control problem can be 
formulated in a condensed form as 

(P) MinF(/), /GJ’ad. (16) 

Referring to [44, Theorem 3.1], we know that the control problem (P) has 
at least one (optimal) solution / for all z/ > 0 . To determine this solution 
numerically, we need the first and second-order derivatives of the objective 
function F. Since the mapping / i-G (iz, i;)^ is twice continuously differentiable, 
so is F : LP{Q) — > M. Its first derivative F'{f) in the direction h G L^iQ) 
can be computed as follows 


T 

F'{f)h = jj dxdt {{<ff + vf)h} , ( 17 ) 


where (pf denotes the first component of the so-called adjoint state 
It solves a linearized FitzHugh-Nagumo system, backwards in time. 
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-dtifif - Aipf + R'{uf )(pf - 7 V’/ = c^{uf - Ud), 

-dt^pf + /3ipf + aipf = c^{vf -Vd), ^ 

with homogeneous Neumann-boundary and terminal conditions (ff{x,T) = 
c^{uf{x, T) — ut{x)) and '0/(x, T) = c^{vf{x, T) — vt{x)) in Q. 

This first derivative is used in numerical methods of gradient type. Higher 
order methods of Newton type need also the second derivative F"{f). It reads 

F"{f)h^ = jdx {4r,h{x,Tf + 4(:h{x,Tf} 

Q 

T T (19) 

+ JJ dxdt {[cd - + + v JJ dxdt {h^} , 

on on 

in a single direction h G Lp{Q). In this expression, the adjoint state {r]h, Ch) • = 
G'{f)h denotes the solution of a linearized FitzHugh-Nagumo system similar 
to Eq. (18), see Ref. [44, Theorem 2.2] for more information. 

2.2 First-Order Optimality Conditions 

We emphasize that the control problem (P) is not necessarily convex. Al¬ 
though the objective function J{u,v, f) is convex, in general the nonlinearity 
of the mapping / i-G (uf^Vf)^ will lead to a non-convex function, F. There¬ 
fore, (P) is a problem of non-convex optimization, possibly leading to several 
local minima instead of a single global minimum. 

As in standard calculus, we invoke first-order necessary optimality con¬ 
ditions to find a (locally) optimal control /, denoted by /. In the case of 
unconstrained control, i.e., F^d •= fhe first derivative of F must be 

zero, F'{f) = 0. Computationally, this condition is better expressed in the 
weak formulation 

T 

F'{f)f = JJ dxdt {(^ + =0, V/ e LP{Q), (20) 

where (f denotes the first component of the adjoint state associated with /. If 
/ is not locally optimal, one finds a descent direction d such that F'{f)d < 0. 
This is used for methods of gradient type. 

If the restrictions F^d are given by Eq. (15), then Eq. (20) does not hold 
true in general. Instead, the variational inequality 

T 

F' if) (/-/) = JJ dxdt {i^ + Uf) (/-/)}> 0, V/ e TVd n Uif) , (21) 

on 

must be fulfilled, cf. [45]. Here U{f) C Lp{Q) denotes a neighborhood of 
/. Roughly speaking, it says that in a local minimum we cannot find an 
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admissible direction of descent. A gradient method would stop in such a point. 
A pointwise discussion of Eq. (21) leads to the following identity: 

f{x,t) = , for zy > 0. (22) 

Here, Proj^j^ (x) = min{max{/a, x},denotes the projection to the in¬ 
terval [/a, A] such that /(x,t) belongs to the set of admissible controls Aid 
defined in Eq. (15). According to Eq. (22), as long as ip does not vanish, a 
decreasing value of z/ > 0 yields an optimal control growing in amplitude until 
it attains its bounds A fb^ respectively. Thus, the variational inequality 
Eq. (21) leads to so-called bang-bang-controls [45] for z/ = 0 and p 0. These 
are control signals which attain its maximally or minimally possible values for 
all times and everywhere in the spatial domain i?. A notable exception is the 
case of exactly realizable desired trajectories and z/ = 0, already discussed in 
Sect. 1. In this case, it can be shown that p vanishes [32] and Eq. (22) cannot 
be used to determine the control signal /. 

Numerically, solutions to optimal control are obtained by solving the con¬ 
trolled RD system Eq. (13) and the adjoint system, Eq. (18), such that the 
last identity, Eq. (22), is fulfilled. In numerical computations with very large 
or even missing bounds A? A? Eq. (22) becomes ill-conditioned if u is close to 
zero. This might lead to large roundoff errors in the computation of the control 
signal and can affect the stability of numerical optimal control algorithms. 


2.3 Example 1: Analytical and Optimal Position Control 

In 1972, Schlogl discussed the auto-catalytic trimolecular RD scheme [1,38] 
as a prototype of a non-equilibrium first order phase transition. The reaction 
kinetics R for the chemical with concentration u{xR) is cubic and can be 
casted into this dimensional form R{u) = u{u — a){u — 1). The associated 
controlled RD equation reads 

dtu — d^u u {u — a) {u — 1) = /(x, t), 0 < a < 1, 

in one spatial dimension, x G M. A linear stability analysis of the uncontrolled 
system reveals that iz = 0 and u = 1 are spatially homogeneous stable steady 
states (HSS) while iz = a is an unstable homogeneous steady state. In an 
infinite one-dimensional domain, the Schlogl model possesses a stable traveling 
front solution whose profile is given by 

C/c(e = l/(l+exp(e/\/2))> (23) 

in the frame of reference ^ = x — ct co-moving with front velocity c. This front 
solution establishes a heteroclinic connection between the two stable HSS for 
^ ^ Too and travels with a velocity c = (1 — 2 a) /^/2 from the left to the 
right. 
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As an example, we aim to accelerate a traveling front according to the 
following protocol of motion 

0(t) = -10 + ct H- 200 ^ 

while keeping the front profile as close as possible to the uncontrolled one. In 
other words, our desired trajectory reads Ud{x^t) = Uc{x — 0(t)), and conse¬ 
quently the initial condition of both the controlled and the desired trajectory 
are uq{x) = U(i{x^ 0) = Uc{x-\-lQi). In our numerical simulations, we set T = 20 
for the terminal time T, i? = (—25, 25) for the spatial domain, and the thresh¬ 
old parameter is kept fixed at a = 9/20. Additionally, we choose the terminal 
state to be equal to the desired trajectory, ut{x) = Ud{x,T), and set the re¬ 
maining weighting coefficients to unity, = 1, in the optimal control 

problem. The space-time plot of the desired trajectory Ud is presented in Fig. 
la for the protocol of motion 0(t) given by Eq. (24). 

Below, we compare the numerically obtained solution to the optimal con¬ 
trol problem (P) with the analytical solution from Sect. 1 for the Schlogl 
model. The Schlogl model arises from Eq. (12) by setting a = 0 and ignor¬ 
ing the inhibitor variable v. Consequently, all weighting coefficients associated 
with the inhibitor trajectory are set to zero, = 0, in the functional 

J, Eq. (14). 

Eig. lb depicts the solution for the analytical position control /an which is 
valid for a vanishing Tikhonov regularization parameter u = 0. The numeri¬ 
cally obtained optimal control f for u = 10“^, shown in Eig. Ic, does not differ 
visually from the analytic one. Both are located at the front position where 



(a) (b) (c) 


Fig. 1. (a) Space-time plot of the desired trajectory Ud{x,t) = Uc{x — /(t)) with 
the protocol of motion 0(t) as given in Eq. (24), (b) analytic position control signal 
fan{x,t), Eq. (12), and (c) numerically obtained optimal control / for Tikhonov 
regularization parameter u = 10“^ are presented. The magnitude of the control 
signal is color-coded. In the center panel (b), the dashed line represents /(t). The 
remaining parameter values are a = 9/20, T = 20, and — 1. 
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the slope is maximal, Ud = 0.5 (dashed line in Fig. lb), and their magnitudes 
grow proportional to (j){t). For a quantitative comparison, we compute the 
distance between analytical and optimal control signal ||/ — /an ||2 in the sense 
of L^((5), Eq. (3), and normalized it by the size of the space-time-cylinder 
\Q\=T\Q\ 

ll^lb := |^||^||L2(g). (25) 

The top row of Tab. 1 displays the distance ||/ — /an ||2 as a function of the 
regularization parameter z/. Even for a large value of z/ = 1, the distance is less 
than 5 x 10“^. Decreasing the value of z/ results in a shrinking distance ||/ — 
/an ||2 until it saturates at 8 x 10“^. The saturation is due to numerical and 
systematic errors. Numerical computations are affected by errors arising in the 
discretization of the spatio-temporal domain and the amplification of roundoff 
errors by the ill-conditioned expression for the control, Eq. (22). A systematic 
error arises because the optimal control is computed for a bounded interval 
i? = (—25,25) with homogeneous Neumann-boundary conditions while the 
analytical result is valid only for an infinite domain. 

Another interesting question is how close the controlled state Uf ap¬ 
proaches the desired trajectory Ud. The bottom row of Tab. 1 shows the 
distance between the optimal controlled state trajectory Uf and the desired 
trajectory for different values of u. Similarly as for the control signal, the 
difference lessens with decreasing values of u. Note that the value does not 
saturate and becomes much smaller than the corresponding value for the dif¬ 
ference between control signals. Here, no discretization errors arise because 
a discretized version of the desired trajectory is used as the target distribu¬ 
tion. Nevertheless, systematic errors arise because neither the initial and final 
desired state nor the desired trajectory obey Neumann-boundary conditions. 
This results in an optimal control signal exhibiting bumps close to the domain 
boundaries. However, the violation of boundary conditions can be reduced by 
specifically designed protocols of motion. The further the protocol of motion 
keeps the controlled front away from any domain boundary the smaller is the 
violation of homogeneous Neumann-boundary conditions since the derivatives 
of traveling front solution Eq. (23) decay exponentially for large \x\. An alter¬ 
native way of rigorously avoiding artifacts due to the violation of boundary 


Table 1. The distance ||/—/an ||2 between the analytical control signal /an, valid for 
= 0, and the optimal control / obtained numerically for finite > 0 decreases with 
decreasing values of z/ (top row). Similarly, the optimally controlled state trajectory 
Uf approaches the desired trajectory measured by \\uf — Ud\\ 2 , for smaller values 
of ly (bottom row). 


ly 

1 

E-1 

E-2 

E-3 

E-4 

E-5 

E-6 

ll/-/an||2 

\\uf -Udh 

4.57E-4 

4.77E-4 

1.14E-4 

7.49E-5 

2.50E-5 

8.34E-6 

l.OlE-5 

8.52E-7 

8.40E-6 

8.55E-8 

8.30E-6 

8.56E-9 

8.29E-6 

8.56E-10 
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conditions is the introduction of additional control terms acting on the domain 
boundaries, see Ref. [31]. 

For the example discussed above the numerical optimal control / for z/ > 0 
is computed with a Newton-Raphson-type root finding algorithm. This itera¬ 
tive algorithm relies on an initial guess for the control signal, which is often 
chosen to be random or uniform in space. The closer the initial guess is to the 
final solution, the fewer steps are necessary for the Newton-Raphson method 
to converge on the final solution. The similarity of the numerical and ana¬ 
lytical control solution, see Fig. 1 and Tab. 1, motivates the utilization of 
the analytical result /an as an initial guess in numerical algorithms. Even for 
a simple single component RD system defined on a relatively small spatio- 
temporal domain Q as discussed in this section, the computational speedup 
is substantial. The algorithm requires only 2/3 of the computation time com¬ 
pared to random or uniform starting values for the control. In particular, we 
expect even larger speedups for simulations with larger domain sizes. 


3 Sparse Optimal Control 

In applications, it might be desirable to have localized controls acting only in 
some sub-areas of the domain. So-called sparse optimal controls provide such 
solutions without any a priori knowledge of these sub-areas. They result in 
a natural way because the control has the best influence in these regions to 
achieve a certain objective functional to be minimized. 

For inverse problems, it has been observed that the use of an L^-term in 
addition to the L^-regularization leads to sparsity [46-48]. The idea to use the 
I/^-term goes back to Ref. [49]. 

To our knowledge, sparse optimal controls were first discussed in the con¬ 
text of optimal control in Ref. [34]. In that paper, an elliptic linear model 
was discussed. Several publications followed, investigating semi-linear elliptic 
equations, parabolic linear, and parabolic semi-linear equations; we refer for 
instance to Refs. [35-37] among others. 

In this section, we follow the lines of Refs. [44, 50] and recall the most 
important results for the sparse optimal control of the Schlogl-model and the 
FitzHugh-Nagumo equation. 

3.1 The Control Problem 

In optimal control, sparsity is obtained by extending the objective functional 
J by a multiple of j{f) := H/HLRg)? Ihe L^-norm of the control /. Therefore, 
recalling that J{uf,Vf,f) =: T{f), we consider the problem 

(Psp) Min^(/) + /.j(/), /G J’ad. 

for /^ > 0. The first part F of the objective functional is differentiable while 
the I/^-part is not. 
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Our goal is not only to derive first-order optimality conditions as in the 
previous section but also to observe the behavior of the optimal solutions for 
increasing k and u is tending to zero. For that task we also need to introduce 
second-order optimality conditions. 

As before, there exists at least one locally optimal solution / to the problem 
(Psp), denoted by /. We refer to Ref. [44, Theorem 3.1] for more details. While 
T is twice continuously differentiable, the second part j(/) is only Lipschitz 
convex but not differentiable. For that reason, we need the so-called subdif¬ 
ferential of j(/). By subdifferential calculus and using directional derivatives 
of j{f), we are able to derive necessary optimality conditions. 


3.2 First-Order Optimality Conditions 

We recall some results from Refs. [44,50]. Due to the presence of j{f) in 
the objective functional, there exists a A G dj{f) such that the variational 
inequality Eq. (21) changes to 

T 

JJ dxdt {{<f + iyf + K\){f - f)} >0, V/G Jad nf/(/). (26) 

0i7 

For the problem (Psp), a detailed and extensive discussion of the first-order 
necessary optimality condition leads to very interesting conclusions, namely 

/(x,t) = 0, if and only if |(^(x,t)| < /^, (27) 

f{x,t) = Proj[^_^ (^-^[^{x,t) + KX{x,t)]^ , (28) 

X{x,t) = Proj[_i + 1 ] , (29) 

if z/ > 0. We refer to Refs. [36, Corollary 3.2] and [51, Theorem 3.1] in which 
the case z/ = 0 is discussed as well. 

The relation in Eq. (27) leads to the sparsity of the (locally) optimal 
solution /, depending on the sparsity parameter k. In particular, the larger 
the choice of k is the smaller does the support of / become. To be more 
precise, there exists a value kq > oo such that for every k > kq the only local 
minimum / is equal to zero. Obviously, this case is ridiculous and thus, one 
needs some intuition to find a suitable value k. We emphasize that A is unique, 
see Eq. (29), which is important for numerical calculations. 


3.3 Example 2: Optimal and Sparse Optimal Position Control 

Eor the numerical computations, we follow the lines of Ref. [44] and use a 
non-linear conjugate gradient method. The advantage of using a (conjugate) 
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gradient method lies in the simplicity in its implementation and in the ro¬ 
bustness of the method to errors in the solution process. Moreover, it allows 
to solve the systems Eq. (13) and the adjoint system separately. The disad¬ 
vantage is clearly the fact that it might cause a huge amount of iterations to 
converge, cf. Ref. [44, section 4]. 

Hence, we modify our approach by the use of Model Predictive Control 
[52,53]. The idea is quite simple: Instead of optimizing the whole time-horizon, 
we only take a very small number of time steps, formulate a sub-problem and 
solve it. Then, the first computed time-step of the solution / of this smaller 
problem is accepted on [0, ti] and is fixed. A new sub-problem is defined by 
going one time-step further and so on. Although the control gained in this 
way is only sub-optimal, it leads to a much better convergence-behavior in 
many computations. 

Next, we revisit the task to extinguish a spiral wave by controlling its tip 
dynamics such that the whole pattern moves out of the spatial domain towards 
the Neumann boundaries [9,21,54]. To this goal, following Ex. 6 from Ref. [44, 
Section 4], we set the protocol of motion to 0(t) = (0, min{120,1/161})^ and 
Ud{x,t) := where i^nat denotes the naturally develop spiral 

wave solution of the activator u to Eq. (13) for / = 0. In our numerical 
simulation, we take only 4 time-steps in each sub-problem of the receding 
horizon. Moreover, we set the kinetic parameters in the EHN model, Eq. (13), 
to a = 0.005, Of = 1, yd = 0.01, 7 = 0.0075, and ^ = 0. Eurther, we fix the 
simulation domain Q = ( — 120,120) x ( — 120,120), the terminal time T = 
2000, n = 10“^ as Tikhonov parameter, and fa = —^ and /^ = 5 as bounds 
for the control, respectively. As initial states {uq^vq)'^ a naturally developed 
spiral wave whose core is located at (0, 0) is used; uq is presented in Eig. 2a. 

In addition, an observation-function G L^{Q) instead of the constant 
factor G M is used with a support restricted to the area close to the desired 
spiral-tip. To be more precise, c^(x,t) = 1 holds only in the area defined by 
all (x,t) G Q such that \x — x{t)\ < 20 and vanishes identically otherwise. 
The other coefficients c^, c^, and are set equal to zero. 

The reason for the choice of such an observation-function is clear: a most 
intriguing property of spiral waves is that despite being propagating waves 
affecting all accessible space, they behave as effectively localized particles-like 
objects [55]. The particle-like behavior of spirals corresponds to an effective 
localization of so called response functions [56,57]. The asymptotic theory of 
the spiral wave drift [58] is based on the idea of summation of elementary 
responses of the spiral wave core position and rotation phase to elementary 
perturbations of different modalities and at different times and places. This 
is mathematically expressed in terms of the response functions. They decay 
quickly with distance from the spiral wave core and are almost equal to zero 
in the region where the spiral wave is insensitive to small perturbations. 

The numerical results for the sup-optimal control {k = 0) and for the 
sparse sub-optimal control (/^ = 1) are depicted in Eig. 2b and Eig. 2c, re- 
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Fig. 2. (a) Spiral wave solution of the activator u to Eq. (13) with / = 0. The 
latter is used as initial state uq{x) = r^nat(cc — </>( 0 ), 0 ) in the problem (Psp)- ( b ) 
Numerically obtained sub-optimal control (/^ = 0) and (c) sparse sub-optimal con¬ 
trol solution — 1 ), both shown in the X 2 -t-plane for = 0 with associated 
spiral-tip trajectory (black line). The magnitude of the control signal is color-coded. 
The remaining system parameters are a — 0.005, (a = l, /3 = 0.01, 7 = 0.0075, and 
(5 = 0. In the optimal control algorithms, we set y — 10~®, fa — —5, and fb = 5. 


spectively. One notices that the prescribed spiral tip trajectory is realized for 
both choices for the sparsity parameter k, viz., k = 0 and k = 1. The traces of 
the spiral tip is indicated by the solid lines in both panels. Since the spiral tip 
rotates rigidly around the spiral core which moves itself on a straight line ac¬ 
cording to 0(t), one observes a periodic motion of the tip in the X 2 -t-plane. 
In addition, the area of non-zero control (colored areas) is obviously much 
smaller for non-zero sparsity parameter k compared to the case k = 0 , cf. 
Fig. 2b and Fig. 2c. However, in this example we observed that the amplitude 
of the sparse control to be twice as large compared to optimal control {hz = 0). 

3.4 Second-Order Optimality Conditions and Numerical Stability 

To avoid this subsection to become too technical, we only state the main re¬ 
sults from Ref. [50]. We know for an unconstrained problem with differentiable 
objective-functional that it is sufficient to show F'{f) = 0 and F"{f) > 0 to 
derive that / is a local minimizer of F, if F were a real-valued function of one 
real variable. More details about the importance of second order optimality 
conditions in the context of PDF control can be found in Ref. [59]. 

In our setting, considering all directions h 0 out of a certain so-called 
critical cone (7/, the condition for u > 0 reads 

F"(/)/i2 > 0. 

Then, / is a locally optimal solution of (Psp)- The detailed structure of Cf is 
described in Ref. [50]; also the much more complicated case = 0 is discussed 
there. 
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The second-order sufficient optimality conditions are the basis for interest¬ 
ing questions, e.g. the stability of solutions for perturbed desired trajectories 
and desired states [50]. Moreover, we study the limiting case of Tikhonov 
parameter v tending to zero. 


3.5 Tikhonov parameter tending to zero 

In this section, we investigate the behavior of a sequence of optimal controls 
and the corresponding states as solutions of the problem (Psp) as z/ ^ 0. 
For this reason, we denote our control problem (Pzy), the associated optimal 
control with and its associated states with {uy, Vy) for a fixed z/ > 0. Since 
is bounded in L^{Q), any sequence of solutions {fy}y>o of (Py) has 
subsequences converging weakly* in L^{Q). For a direct numerical approach, 
this is useless but we can deduce interesting consequences of this convergence 
using second order sufficient optimality conditions. 

Assume that the second order sufficient optimality conditions of Ref. [50, 
Theorem 4.7] are satisfied. Then, we derive a Holder rate of convergence for 
the states ^ 

lim —j= {||u^ - wo||l2(q) + Wvu - vo||l2(q)} = 0, (30) 

>-U yj V 

with {uy^Vy) = G{fy) aud (^ 20 , 1 ) 0 ) = G(/o)- We should mention that this 
estimate is fairly pessimistic. All of our numerical tests show that the con¬ 
vergence rate is of order z/, i.e., we observe a Lipschitz rather than a Holder 
estimate [50]. As mentioned in Ref. [50], it should also be possible to prove 
Lipschitz stability and hence to confirm the linear rate of convergence for 
z/ ^ 0 with a remarkable amount of effort. 


3.6 Example 3: Sparse Optimal Control with Tikhonov parameter 
tending to zero 

Finally, we consider a traveling pulse solution in the FitzHugh-Nagumo system 
in one spatial dimension N = 1. Here, the limiting case of vanishing Tikhonov 
parameter, z/ = 0, is of our special interest. We observe that Newton-type 
methods yield very high accuracy even for very small values of z/ > 0. This 
allows us to study the convergence behavior of solutions for z/ tending to zero 
as well. 

Following Ref. [50] and in contrast to the last example in Sect. 3.3, we solve 
the full forward-backward-system of optimality. We stress that this is numer¬ 
ically possible solely for non-vanishing value of z/. However, we constructed 
examples where an exact solution of the optimality system for z/ = 0 is acces¬ 
sible as shown in Ref. [50, Section 5.3]. In this sequel, our reference-solution, 
denoted by tZref, will be the solution of (P^y) for z/ := z/ref = 10“^^. For smaller 
values the numerical errors do not allow to observe a further convergence. The 
distance \\uy — Uref\\L‘^(Q) stagnates between u = 10“^^ and z/ref < 10“^^. 
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Next, we treat the well-studied problem of pulse nucleation [60, 61] by 
sparse optimal control. We aim to start and to stay in the lower HSS for the 
first two time-units, i.e., Ud{x,t) = —1.3 for t G (0,2). Then, the activator 
state shall coincide instantaneously with the traveling pulse solution he., 
= Unat{x,t — 2). To get the activator profile we solve Eq. (13) 
for / = 0 and its profile is shown in Fig. 3a. 

In our optimal control algorithms, we set the parameters to i? = (0,75), 
T = 10, (a = 1, /3 = 0, 7 = 0.33, and 5 = —0.429. Moreover, here we use 
a slightly different nonlinear reaction kinetics R{u) = u{u — \/3){u + \/3) in 
Eq. (13) but this does not change the analytical results. The upper and lower 
bounds for the control / are set to very large values, viz. fa = —100 and 
fb = 100. In addition, the coefficients in Eq. (14) are kept fixed, viz., = 1 
and = 0. 



Fig. 3. (a) Segment of a traveling pulse solution (i^o,r'o)^ in the uncontrolled 
FitzHugh-Nagumo system, Eq. (13) with / = 0. (b) Numerically obtained sparse 
optimal control solution for almost vanishing Tikhonov parameter, u = 10“^°, 
and (c) the associated optimal state Uj^. The amplitude of the control signal is color- 
coded. The kinetic parameter values are set too = 1, ^5 = 0,7 = 0.33, (5 = —0.429, 
and R{u) = u{u — ^){u + a/3). 


Our numerical results obtained for a sparse optimal control acting solely 
on the activator cf. Eq. (13), are presented in Fig. 3b and Fig. 3c. In order 
to create a traveling pulse solution from the HSS Ud = —1.3, the optimal 
control resembles a step-like excitation with high amplitude at x 40. Since 
the Tikhonov parameter is set to = 10“^^, large control amplitudes are to be 
expected and indicate that in the unregularized case, even a delta distribution 
might appear. Because this excitation is supercritical a new pulse will nucleate. 
In order to inhibit the propagation of this nucleated pulse to the left, the 
control must act at the back of the pulse as well. Thus, we observes a negative 
control amplitude acting in the back of the traveling pulse. We emphasize 
that the desired shape of the pulse is achieved qualitatively. The realization 
of the exact desired profile can not be expected due to a non-vanishing sparse 
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parameter k, = 0.01. Even for this respectively small value, the sparsity of the 
optimal control shows. 

Since the displayed control and state are computed for an almost vanishing 
value z/ = 10“^^, we take the associated state as reference-state Uref in order 
to study the dependence of the distance \\uu — '^ref||L 2 (g) on z/ > 0. From 
Table 2, one notices the already mentioned Lipschitz-rate of convergence for 
decreasing values z/ > 0, \\uu — '^ref||L2(g) oc z/. This observation is consistent 
to results from [50] for various other examples. 


Table 2. The comparison of the distance \\uu — tZref||L2(g) between the numerically 
obtained states and the numerically obtained reference-solution tzref, computed 
for ly = 10“^°, for decreasing values of z/ > 0. 


ly 

lE-3 

lE-4 

lE-5 

lE-6 

lE-7 

lE-8 

lE-9 

\\Uu — tiref |L2(g) 

1.58E-1 

1.16E-2 

1.33E-3 

1.35E-4 

1.35E-5 

1.34E-6 

1.31E-7 


4 Conclusion 

Optimal control of traveling wave patterns in RD systems according to a 
prescribed desired distribution is important for many applications. 

Analytical solutions to an unregularized optimal control problem can be 
obtained with ease from the approach presented in Sect. 1. In particular, the 
control signal can be obtained without full knowledge about the underlying 
nonlinear reaction kinetics in case of position control. Moreover, they are a 
good initial guess for the numerical solution of regularized optimal control 
problems with small regularization parameter z/ > 0, thereby achieving a 
substantial computational speedup as discussed in Sect. 2.3. Generally, the 
analytical expressions may serve as consistency checks for numerical optimal 
control algorithms. 

For the position control of fronts, pulses, and spiral waves, the control 
signal is spatially localized. By applying sparse optimal position control to 
react ion-diffusion systems, as discussed in Sect. 3, the size of the domains 
with non-vanishing control signals can be further decreased. Importantly, the 
method determines sparse controls without any a priori knowledge about re¬ 
strictions to certain subdomains. Additionally, sparse control allows to study 
second order optimality conditions that are not only interesting from the the¬ 
oretical perspective but also for numerical Newton-type algorithms. 
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